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1 Introduction 

Pair production processes for electroweak vector bosons provide a rich spectrum of observ¬ 
ables, which are crucial to test in depth the SU{2)l x U{1)y gange sector of the Standard 
Model. In particular, the production of pairs of resonant vector bosons allows for precise 
studies of the electroweak triple gauge couplings, while considering off-shell vector boson 
pairs is required for precision Higgs phenomenology. Fnrthermore, diboson production 
processes are important backgronnds in direct new physics searches. The main production 
channel for pairs of vector bosons at hadron colliders is quark-antiquark annihilation and 
great progress has been achieved in the last years with the compntation of the next-to-next- 
to-leading order (NNLO) QCD corrections to qq —)• 77 [1], qq —)• Z 7 [2], qq —)• ZZ [3] and 
qq —>■ W^W~ [4] production at the LHC. Furthermore, the fermionic NNLO corrections 
to qq —)• 7 * 7 * were derived in [5]. 

The gluon fusion channel contribntes to 77 , ZZ, Z 7 and W~^W~ production. As 
a quark-loop induced process, its leading order (LO) cross section is suppressed by two 
powers of with respect to that of the quark channel. This implies that it formally 
contribntes only at NNLO in the perturbative expansion of the hadronic process, but 
numerical enhancements may be expected due to the large gluon luminosities at typical 
energies for diboson production at the LHC. For the gluon-induced processes, the one- 
loop amplitudes and the corresponding one-loop squared interference terms have been 
computed long ago [6-12]. Their impact on the total cross section was found to range 
approximately from 5% to more than 10% for different hnal states at the LHC, and to 
rise with increasing collider energy [1-4]. These values can substantially increase up to 
about 30% when particular sets of cuts, relevant for example for Higgs boson searches, 
are applied [13, 14]. It is therefore clear that the inclusion of gluon channel contribntions 
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can be important in order to achieve a description of the full process which matches the 
experimental precision. Beyond the actual size of the known leading order corrections in 
the gluon channel, it is unclear how large the associated theory uncertainty actually is. 
By comparison with Higgs production in gluon fusion [15-17], the conventional LO scale 
variation is not expected to allow for a reliable estimate of the size of neglected higher 
order corrections. The recent NNLO predictions for the total ZZ and W^W~ production 
cross sections take into account the quark channel at NNLO and the gluon channel at 
LO, resulting in a scale uncertainty of about 3% [3, 4]. In order to thoroughly control the 
theory uncertainty to this level of precision, it is therefore very desirable to compute the 
next-to-leading order (NLO) contributions for the gluon induced subprocess. Currently 
this has been done only for gg —>■ 77 [18, 19], and the NLO corrections have been found 
to be not only sizeable but also important for stabilising the theoretical predictions [19]. 
Finally, precise theoretical predictions for gg ^ ZZ can be useful for constraining the total 
Higgs boson decay width at the LHC [20-23]. 

Technically, the computation of the NLO corrections to gg —)• V 1 V 2 requires two in¬ 
gredients, the two-loop virtual corrections to gg —)• V 1 V 2 and the one-loop real-virtual 
corrections to the corresponding radiative processes with one more parton in the final 
state. By now the computation of the one-loop amplitudes with an extra gluon does not 
constitute any conceptual difficulty and can be pursued with standard techniques for one- 
loop multi-legs processes [24-30]. The two-loop amplitudes, on the other hand, are known 
only for gg —?■ 77 [18] and for gg —?■ Z 7 [31], in both cases for on-shell final state photons. 
In order to obtain physical predictions, both contributions need to be combined using a 
subtraction scheme to isolate and cancel unphysical IR divergences. In this case, a NLO 
scheme [32, 33] would be sufficient. 

In this paper we calculate the missing two-loop massless QCD corrections to gg —)• 
V 1 V 2 , with an off-shell vector boson pair V 1 V 2 = , ZZ, Z^*,W^W~. The calculation 

builds upon the master integrals for four-point functions with massless propagators and two 
massive external legs, which were computed recently in the case of equal masses in [34, 35], 
and in the case of different masses in [36-39]. The former were used for the hrst NNLO 
fully-inclusive calculations oi ZZ [3] and W~^W~ [4] production at the LHC, while the 
latter allowed the computation of the two-loop corrections to —)• V 1 V 2 [39, 40]. A subset 
of these master integrals was also computed independently in [5, 41]. While the inclusion 
of massive top-loop mediated subprocesses would be of interest for some phenomenological 
applications [20, 42], the computation of the two-loop amplitudes requires knowledge of 
challenging new master integrals, which should be addressed in the future. 

The paper is structured as follows. In Section 2 we describe the tensor decomposition 
of the partonic current for the process gg ^ V 1 V 2 and consider the possible electroweak 
coupling structures. We include the vector boson decays and describe the helicity ampli¬ 
tudes for the process gg ^ V 1 V 2 ^ ^ leptons in terms of scalar form factors in Section 3. 
The actual calculation of the loop contributions to these form factors is described in Sec¬ 
tion 4, which includes a dicussion of UV renormalisation, IR subtraction and various checks 
we performed on our results. In Section 5 we present numerical results obtained with our 
C++ implementation. Finally, we conclude in Section 6 . In Appendix A we give explicit 
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formulae for obtaining the physical form factors appearing in the helicity amplitudes from 
the original tensor coefficients computed in this paper. We provide computer readable files 
for our analytical results and our C++ code for the numerical evaluation of the amplitudes 
on our VVamp project page on HepForge at http://vvamp.hepforge.org. 

2 Partonic current for gg ^ V 1 V 2 

We consider the production of two massive off-shell vector bosons, V 1 V 2 , in the gluon fusion 
channel, 

9{Pi) + g{P2) —^ Pi(P3) + ^ 2 (^ 4 ), (2.1) 

where V 1 V 2 = 7 * 7 *, ZZ^ Z^*, W^W~. The final states and W^Z instead are 

forbidden by charge conservation. Since the two vector bosons are off-shell we have in the 
general case 

pI=pI = 0, pi> 0 , P 4 > 0 , pI^pI, ( 2 . 2 ) 

with the usual Mandelstam invariants defined as 

s = {pi+P2f, t = {pi-pzf, u = {p2-pzf, (2.3) 

and the relation 

s -\-1 -\- u = -\- P 4 . (2-4) 

The physical region for the scattering kinematics has the boundary tu = p\p\ and fulfils 

s > {\flii + , \{pI+ pI - s - k) <t <]^{pI+ pI - s + k) (2.5) 

where k is the Kallen function 

k{s,pI,pI) = +p^ + pj-2{spl+plpl+pls). ( 2 . 6 ) 

We denote the scattering amplitude for the process (2.1) by 

S{pi,P2,P3) = Sf,^pa {pi,P2,Ps)ei {pi ) €2 {p2 ) 63 ^ {p3 ) 64 (p4 ) 

where ei, 62 are the polarisation vectors of the incoming gluons, 63 , €4 are the polarisation 
vectors of the outgoing massive vector bosons. Pa = pi + P 2 — P 3 and an overall factor 
is kept implicit with e being the positron charge. Since we will consider leptonic decays 
of the massive vector bosons we will be able to construct the full amplitude including the 
decays from the partonic current 

Spu {pi ,P2,P3) = Spy pa (pi, P 2 , P3 ) e? {Pl ) ^2 {P2 ) 

for the 2 —)• 2 process. In particular, it is only the latter which receives (pure) QCD 
corrections at any order in perturbation theory. 
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In order to compute the partonic current it is useful to consider its tensor decomposi¬ 
tion. Based on Lorentz invariance only, there are 138 independent tensor structures which 
can contribute 




+ V f 6 p^- p^- + 0^^ p^- p 


iij 2 =i 


^ 31J2 y Oi ^^32 ^ ^3132 y i^3l ^32 ^ ^3132 ^ Oi ^^32) 


3 

+ Z] ^hhnhPhPhPhP'^H^ ( 2 . 7 ) 

3lp2 J3 J4 = l 

where the coefficients aj, b^- and Cijki are scalar functions of the kinematic invariants s, t, 
p|, p1 and of the space-time dimension d. Not all structures are relevant for our calculation. 
Many of them simply drop due to the transversality of the gluons’ polarisation vectors 

ei - pi = 62 - ^2 = 0 . ( 2 . 8 ) 

Moreover the tensor structure can be further simplified by fixing explicitly the gauge for 
the incoming gluons. A particularly simple choice is given by the symmetrical condition 


ei • J>2 = 62 • Pi = 0 , 


(2.9) 


which corresponds to the following rules for the polarisation sums 

^ Pl‘P2 

+ (2.10) 

Further conditions can be applied on the polarisation vectors of the massive vector bosons 
hi, V 2 . We employ for their polarisation vectors 


63 • P3 = 64 • P4 = 0 . 


( 2 . 11 ) 


and for the polarisation sums 


ECA3fe)'SA,fe) = -9'“' + 


Aa 


E'rA.(Pi)eSA,(P4) = -9'“' + 

A 4 


M V 

pm 
2 ’ 
pi 

PaPI 

pI 


( 2 . 12 ) 


Imposing the constraints (2.8), (2.9) and (2.11) one is left with only 20 independent 
tensor structures and we can write the partonic current according to 

20 

S^''{pi,P2,P3) = ^ Aj{s,t,pl,pl)Tf'' , (2.13) 


i=i 
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where the Aj are scalar functions of s, t, p\ and d. The tensors are defined as 
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(2.14) 


We stress that the tensor decomposition (2.13) is based only on Lorentz symmetry, gauge 
invariance and the properties of the boson decays and holds therefore at every order in 
perturbative QCD. Moreover, no assumption has been made on the dimensionality of space- 
time and the result is valid for any values of the parameter d. 

The scalar form factors Aj can be extracted from the amplitude (2.13) by applying 
suitable projecting operators. The projectors themselves can be decomposed in the same 
20 tensors as 

20 

for J = 1,..., 20, (2.15) 

i=l 

where also Bji are functions of the external invariants and d. Their explicit form can be 
determined imposing 

P/"' for j = 1,..., 20, (2.16) 

pol 

where the polarisation sums are evaluated in d dimensions according to (2.10) and (2.12). 
The explicit results for the coefficients Bji are rather lengthy and we prefer not to write 
them here explicitly. Computer readable files for the latter are given on our project page 
at HepForge. 

The partonic current is the only one which receives contributions from QCD radiative 
corrections and, for two gluons of helicities Ai and A 2 , can be written as 

Siiu{Pi\P2\P3) = J2 ^v!v 2 SillpaiPl,P2,P3yiXiiPi)^2\2iP2) > ( 2 - 17 ) 

j 

where jg the overall colour structure and the index j runs over different possible 

classes of diagrams discussed below, see also Fig. 1, which are characterised by different 

f?l 

electroweak couplings 

Before proceeding, it is convenient to introduce some notations needed in the following. 
As long as we work in QCD, we only need to consider the coupling of electroweak vector 
bosons V to fermions. We follow [43] and parametrise the couplings as 

i e, where e = -v/dvra is the positron charge , (2-18) 
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the vector bosons couple to the same fermion loop, [A], to different fermion loops, [i?], or to an 
intermediate vector boson, [Fy]- The sum of all type [B] contribution and the sum of all type [Ty] 
contributions vanish, respectively. 


such that all fermion charges are expressed in units of e and 


with 


rt"/i/2 = jV 


/1/2 


1 - 75 


^/l/2 - -^h <5/1/2 


^hf2 - 
rW _ 

-L/f. — 


- sin^6>^,e/^ 
sin cos 6 w 

1 


V1/2 




+ ^/l/2 > ( 

^l + 75\ 

(2.19) 

V 2 j ’ 

II 

1 

e/l <5 /i /2 : 

(2.20) 

II 

1 

sin 0^6/^ 

COsOw ’ 

(2.21) 

Rhh = 0 


(2.22) 


where 6/^/2 is unity for /i / /2, but belonging to the same isospin doublet, and zero 
otherwise. 

Let us consider the different electroweak coupling structures in detail. It is clear that, 
since we do not take any electroweak radiative corrections into account, at least one of 
the two vector bosons must be coupled to an internal fermion loop. In order to compute 
the one- and two-loop massless QCD corrections we need to consider the following three 
possibilities, see Fig. 1. 


Class A: Both vector bosons V1V2 are attached to the same fermion loop. In this case the 

diagrams are proportional to the charge weighted sum of the quark flavours, which 
f^l 

we denote as = Nv-^V2- These diagrams could in principle yield two different 

contributions. One, proportional to the sum of the vector-vector and the axial-axial 
couplings, in which all dependence on 75 cancels out. The second, instead, contains 
the vector-axial coupling and is linear in 75. Due to charge parity conservation this 
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last contribution is expected to always vanish identically for massless quarks running 
in the loops, for any choice of Vi and V 2 [8, 9, 42], One then easily finds that 


i 


i 

^ww = \Y1 


where the indices i,j run over the flavours of the quarks in the loop and Lg.q. = Rq^qi 
such that Nyry = e^.. 

Class B: The two vector bosons are attached to two different fermion loops. This con¬ 
figuration is of course possible only starting from two loops on. Each fermion loop 
contains both a vector and an axial piece. For the case of two-loop massless QCD 
corrections relevant here, both contributions can be shown to vanish. The axial con¬ 
tribution cancels out for degenerate isospin doublets, while the vector piece must sum 
up to zero due to Furry’s theorem. 


Classes Fv: Only for the case of F 1 V 2 = W^W~, one should also take into account the 
s-channel production diagrams, where the incoming gluons produce an intermediate 
electroweak gauge boson V = jZ*, which then decays into the outgoing VF-pair, see 
Fig. 1. Charge-parity invariance ensures that the vector part of these diagrams must 
sum up to zero. Again, the axial part cancels out for degenerate isospin doublets, 
and therefore also in the case of massless quarks running in the loops. 


For the case of the one- and two-loop contributions considered here, we can therefore 
simplify (2.17) to 

St,u{Pi\P2\P3) = NvrV2 Sl^lApi,P2,P3)eiXiiP^)^2X2iP2) > (2-24) 

with Nv^V 2 given in (2.23) and consider the coefficients defined by 

^j{s,t,pl,pl) = 6a,a2NviV2Af\s,t,pl,pl)- (2.25) 

It is instructive to study the transformations of the partonic current (2.24) under 
permutations of the external legs. We define the following two permutations 


7ri2 := Pi ^ P2 ^ {t ^ u}, 

-^34 ■= P3 ^ P4 ^ {t U , Ps^ pI}- (2-26) 


Because of Bose symmetry these two permutations must leave the partonic amplitude 
unchanged. This enforces a well defined behaviour of the coefficients Aj(s, under 

the action of 7ri2 and 7r34. From direct inspection of (2.13) one finds that the following 


- 7 - 




relations must be fulfilled: 

7ri2 : A^f\s,u,pl,pl) = A^f\s,t,pl,pl), 
A^f\s,u,pl,pl) = Alf\s,t,pl,pl), 
A^^\s,u,pI,pI) = A^^\s,t,pl,pl), 
= ^^^\s,t,pl,pl), 
^wis,u,plpl) = A^^\s,t,pl,pl), 
^m\^’'>^^P3iPl) = ^^^\s,t,pl,pl), 
TT34 : a\^\s,u,pI,pI) = ^^^\s,t,pl,pl), 
A^f\s,u,pl,pl) = ^^f\s,t,pl,pl), 
aIj^\s,u,pI,pI) = Alf\s,t,pl,pl), 
A^^\s,u,pI,pI) = -^Sf (s,i,P3iP4) 
A^^\s,u,pI,pI) = -A^£{s,t,phpl) 
A^^\s,u,pI,pI) = ^^^\s,t,pl,pl), 
aI^\s,u,pI,pI) = A\^\s,t,pl,pl), 


a\^\s,u,pI,pI) = Alf^{s,t,pl,pl) , 
A^^\s,u,pI,pI) = A^^\s,t,pl,pl) , 
A^^\s,u,pI,pI) = A[^\s,t,pl,pl) , 
A^^\s,u,pI,pI) = A[^\s,t,pl,pl) , 
A^f^{s,u,pl,pl) = A\^^{s,t,pl,pl) , 

(2.27) 

A\f'{s,u,pl,pl) = ^^^\s,t,pl,pl) , 
A^f'{s,u,pl,pl) = ^^f'{s,t,phpl) , 
A^f'{s,u,pl,pl) = > 

^^S{s,u,pI,pI) = -A^^}{sA,pI,pI) , 

^^iq{s^u,pI,pI) = A^^^{s,t,pl,pl) , 
^is{s,u,pI,pI) = J^^{s,t,pl,pl) , 


(2.28) 

It is interesting to notice that, upon exploiting all of these crossing relations, only 9 out of 
the 20 coefficients Al^^ turn out to be effectively independent, while the other 11 coefficients 
can be obtained by crossing of the external legs. 


3 Helicity amplitudes for gg ^ V 1 V 2 ^ A leptons 


We consider physical processes, where the two off-shell vector bosons decay into lepton 
pairs 

5 ( 7 * 1 ) + 5 ( 7 * 2 ) W(P3) + ^ 2 (^ 4 ) h{Pb) + k{p&) + hip?) + kiPs) (3.1) 

such that p3 = P5 + P 6 , Pa = Pj + Ps and p\ = p^ = Pj = p\ = 0. As long as we 

consider QCD radiative corrections the amplitudes can be written, at any order 

in perturbation theory, as the product of the partonic current for gg —)■ V 1 V 2 with the 
two leptonic currents for the decay products, Vi —)• kk and V 2 —^ hk-, mediated by the 

propagators of the two off-shell vector bosons Vi and V 2 - We write the propagator for an 

off-shell vector boson in the gauge as 


with 



(3.2) 


(3.3) 


D^* {q) = , 


Dz,w{q) = ( 5 ^ - my + iTymy) 


(3.4) 



where my is its mass and Fy is its decay width. In our case the massive vector bosons 
couple to massless fermion lines such that the term proportional to (1 — 0 can be dropped. 

In the following we consider fixed helicities of the external particles and compute the 
amplitudes for the different helicity configurations. While the quantities and formulae 
presented up to this point were treated in d dimensions throughout, we now consider 4- 
dimensional external states in order to compute the amplitudes for specihc helicities. Since 
the decay leptons are massless, helicity is conserved along the leptonic decay currents and 
the amplitude can be written as 

-^AiA2A3A4(Fl>F2;P5,P6,P7,P8), (3.5) 

where Ai and A 2 are the helicities of the incoming gluons, while A 3 and A 4 are the helicities 
of the two leptonic currents. It is clear that there are 16 different helicity conhgurations, 
depending on the different possibilities for the initial and final states. Each gluon has two 
possible helicity states, which we denote by L (-) and R (+), and similarly each leptonic 
current occurs in either left- or right-handed configuration, again denoted by L and i?, 
respectively, such that \j = L,R, for j = 1,...,4. As we will show explicitly later on, 
all 16 helicity configurations can be obtained from only two independent ones, by simple 
permutations of the external legs and complex conjugation. We choose as independent 
conhgurations the following two 


-^LLLl(F1>F2;P5,P6,P7,P8) , :P%) ■ 


Viy2 


(3.6) 


With the notations introduced above we write the two independent helicity ampli¬ 
tudes (3.6), up to two loops, as: 


-^AlA2LL(FyF2;P5,P6,P7,P8) = (4^)^ MaiA2Ll(F 1 , P2; P5, P6, P7, Fs) 


tVi jV2 


(3.7) 


where the basic amplitudes Mx^x 2 Ll(pi,P 2 ',P 5 ,Pe,Pr,Ps) are constructed from the partonic 
current (2.17) and the leptonic currents (3.9) according to 

MAi A 2 LL (py F2; P5, P6, P7, P8) = (FI ) ^ 2 X 2 (P ^) (pi, P2, Ps) P'L (Ps ’ Pd ) (P? ’Ps)' 

(3.8) 

The leptonic decay currents do not receive any QCD corrections and are simple tree-level 
objects. They can be easily expressed in the usual spinor-helicity notation [44, 45] as 

L^iPs^Pe) = U-(P5)7''y+(P6) = [6|7''|5) = (5|7^|6], (3.9) 

Lnipt^Pe) =u+{P5)7^V-{pq) = [5|7^|6) = (T^(P5,P6))* =L^(p 6 ,F^)- (3-10) 


Note that, in this case, a permutation of the external momenta is equivalent to a complex 
conjugation of the current and it corresponds to a hip of the helicity L ^ R. 
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Once the tensor decomposition of the partonic current is hxed, it is straight-forward to 
express the two basic helicity amplitudes M-llll and M-lrll in (3-8) in the usual spinor- 
helicity notation [44, 45]. We replace the gluon polarisation vectors according to 




[ 217 ^ 1 ) 
V2[12] ’ 


eiR(Pi) 


( 217 ^ 1 ] 
V2{21) ’ 




[1|712) 

V2[21] ’ 


''2R 


{P2) 


( 1 | 712 ] 

V2(12) ’ 

(3.11) 


which is of course compatible with the polarisation sums (2.10) and (2.12). Note again that 
here we are assuming 4-dimensional external states. This allows to reduce considerably the 
number of independent structures that are required for parametrising a specific helicity 
conhguration. Using (3.11) we find that both basic amplitudes can be written in terms of 
9 independent spinor structures as 

MaiA2Ll(P1,P2;P5,P6,P7,P8) = CaiAz | [2|^3 1) (57) [68 ] 

+ (15)(17)[16] [18] +.^3(15) (27) [16] [28] 

-F (25) (17) [26] [18] -F Sg (25) (27) [26] [28] ) 

+ (15)(17)[16] [28] +^7(15)(17)[26] [18] 

+ eI^^^ (15)(27)[26] [28] -FSg (25)(17)[26] [28] | , (3.12) 

where the 18 newly introduced form factors are simple linear combinations of the 

scalar coefficients Aj. The spinor structure of the amplitudes for the configurations LLLL 
and LRLL differs only by an overall factor which reads in the two cases 

Cll = [1^32)^^, Clr = [2]^3l), (3.13) 

but the form factors Ej"^ and E^^ are different. We also note, in passing, that the spinor 
structure of (3.12) exhibits also a formal similarity to that of the RLL amplitude for 
qq' —)• U 1 V 2 —^ hlehh [39, 40], again up to an overall factor and with, of course, completely 
unrelated form factors. Similar as before, we also define the functions 

E}^^^is,t,plpl) = Sa,a,Nv,v,Ej^^^^^\s,t,plpl). (3.14) 

The explicit expressions for the form factors Ej^^^ in terms of the coefficients Aj are given 
in Appendix A. 

In order to obtain all 16 helicity amplitudes from (3.12), one should recall that complex 
conjugation has the effect of reversing the helicity of the external gluons, 


(eiL(Pi))* = eii?(Pi), i^LiP^))* = 4RiP2) > (3.15) 

and similarly for the leptonic currents, see (3.10) (3.9). We define with the symbol [...]‘^ a 
complex-conjugation operation which, when applied on the amplitudes MaiAzLLj ^-cts only 
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on the spinor structures, i.e. leaves invariant the form factors Given the explicit 

form of (3.12), it is easy to see that this corresponds to simply exchanging angle brackets 
with squared bracket and vice versa 

[MaiAzLl]^ = MaiAsLL ((g) ^ [g]) • (3.16) 

Hence, we can derive the missing helicity amplitudes for left-handed leptonic currents from 
the two basic amplitudes as 

MRLLL{Pl,P2;P5,Pe,P7,P8) = [^LRLl{P1,P2;P6,P5,P8,P7)]^' , 
M.RRLL{Pl,P2]P5,Pe,P7,P8) = [M.LLLL{Pl,P2]P6,P5,P8,P7)f , (3.17) 

where one should note that the lepton and anti-lepton momenta are exchanged in the r.h.s. 
in order to have a left-handed leptonic currents on the l.h.s. The corresponding formulae 
for the basic amplitudes for right-handed leptonic currents can be obtained from the ones 
above by simple permutations of the lepton and anti-lepton momenta 

MaiA2Rl(P1,P2;P5,P6,P7,P8) = MaiA2Ll(P 1 , P2;P6,P5,P7, Ps) , 
MaiA2Lr(P1,P2;P5,P6,P7,1>8) = Mx^^X^lUpi, P2 ; P5, P6, P8, P7) , 
'^XiX2Rr{P 1^P2-,P5,P6,P7,P8) = M.XiX2LL{Pl,P2-,Pe,P5,P8,P7) ■ (3-18) 

With these formulae also all the 16 physical amplitudes in (3.7) can be easily 

obtained, recalling that in the case of right-handed leptonic currents one should, of course, 
exchange the corresponding couplings L^.f. -H- R^^f - 

As we already stated above, the partonic current receives contributions form QCD 
radiative corrections and it can be expanded as 

Sf,upaiPl,P2,P3) = (^) Sj}Jp„{pi,p2,P3) + (^) Sj^Jp„{pi, P2, P 3 ) + 0{al) , (3.19) 

where obviously the perturbative expansion starts only at one-loop order. Of course also 
the coefficients Aj, and equivalently the form factors have the same expansion 



4 Calculation of the form factors 

The calculation of the coefficients Ej^^^ proceeds as follows. We produce all one- and 
two-loop Feynman diagrams relevant for gg —?■ V 1 V 2 using Qgraf [46]. In particular we 
focus only on diagrams in classes A and B with massless quarks, for which we hnd 8 di¬ 
agrams at one loop and 138 diagrams at two loops. Diagrams in class Ey, in fact, are 
simple three-point functions, which sum up to zero due to charge-parity invariance. The 
coefficients Aj are then calculated by applying the projectors defined in (2.15) on the dif¬ 
ferent Feynman diagrams. We insert the Feynman rules in our diagrams, where we employ 
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the Feynman-’t Hooft gauge = 1) for internal gluons. After evaluation of Dirac traces 
and contraction of Lorentz indices every Feynman diagram is expressed as linear combina¬ 
tion of a large number of scalar integrals. The latter belong to the family of the massless 
four-point functions with two off-shell legs of different virtualities and can be reduced to 
a small set of master integrals using integration-by-parts identities [47-50]. We employ 
Reduze 2 [51-54] to map all scalar integrals to the three integral families given in [35] and 
their crossed versions, and subsequently to reduce them to master integrals. In this way, 
we obtain analytical expressions for the coefficients Aj as linear combinations of the lat¬ 
ter. For the master integrals we employ the solutions presented in [39]. With the explicit 
expressions for the coefficients Aj at the different perturbative orders, it is easy to obtain 
the corresponding results for the using the formulae given in Appendix A. Form [55] 

was used extensively for all intermediate algebraic manipulations. 


Because of the lack of any tree-level contribution to the process gg —)• V 1 V 2 , the UV and 
IR pole-structure of the one- and two-loop amplitudes is very simple. Clearly, the one-loop 
amplitude must be both UV- and IR-finite, and therefore the pole structure of the two-loop 
amplitude will be, effectively, what one usually encounters for a one-loop QCD amplitude. 
As discussed above, QCD radiative corrections affect only the partonic amplitude and can 
be taken into account via the 20 independent scalar coefficients Aj, see Eq. (2.13). Working 
in conventional dimensional regularisation, we may alternatively consider the 18 physically 
relevant form factors defined as d dimensional linear combinations of the Aj, see 

Appendix A. In what follows, all considerations regarding UV-renormalisation and the 
structure of the IR poles of the partonic amplitude hold identically for any Aj and 
We will therefore focus on the scalar coefficients rather than on the full partonic amplitude, 
and use the symbol D to refer to any of the latter, 

D G |Aj , , for any j = l,...,20, i = l,...,9, X1X2 = LL,LR. 


We start by performing UV-renormalisation in the MS scheme. In massless QCD this 
amounts to replacing the bare coupling, a^, with the renormalised one, ag = as{E^), where 
g is the renormalisation scale. Here we only need the one-loop relation 


2e o 2e 

ao /^o ^e = Oiag 




(4.1) 


where 


Sg = ( 47 r)^e , with the Euler-Mascheroni constant 7 = 0.5772..., (4.2) 


e = {A—d)/2, go is the mass-parameter introduced in dimensional regularisation to maintain 
a dimensionless coupling in the bare QCD Lagrangian density, and finally /3o is the first 
order of the QCD /3-function 


/3o 


11 Ca — AEp Nj 
6 


with Ca = N , Cp 


iV2 - 1 ^ _ 1 

2N ’ ^ ~ 2' 


(4.3) 
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The renormalisation is performed at = s, the invariant mass squared of the vector-boson 
pair. The renormalised form factors read then, in terms of the un-renormalised ones, 

0(2) = . (4.4) 

After UV renormalisation, the two-loop coefficients O(^) contain still residual IR singu¬ 
larities. In any IR-safe observable these divergences are cancelled by the corresponding ones 
produced in one-loop radiative processes with one more external parton. In the present case 
of gg —)• I4 iV2 , as discussed already above, the IR-poles at two loops are of NLO type and 
their structure has been know for a long time. Here, we choose to follow the conventions 
used for the NNLO corrections to qq ^ V 1 V 2 in [39], which required a NNLO subtraction 
scheme. The exact structure of the IR poles up to NNLO in QCD was predicted first by 
Catani [56]. We present our results in a slightly modified scheme described in [57], which 
is well suited for the gr-subtraction formalism. 

We define the IR finite amplitudes at renormalisation scale g, in terms of the UV 
renormalised ones as follows 


q( 1),finite 
q{ 2),finite 


hW, 


(4.5) 


where for the gluon-fusion channel we have 


/i(e) = /r^*(e) + /r^^(e), 


jsoft 






+ 5 ( 0 ) 

^ qr 


Ca, 


(4.6) 

(4.7) 

(4.8) 


Following [57] we then put 6 q^ = 0. We provide the explicit analytical results for the finite 
remainders of the coefficients Aj in this scheme, obtained for = s, on our project page 
at HepForge. 

Finally, it is straight-forward to convert these finite remainders into the Catani’s origi¬ 
nal subtraction scheme [56] , as extensively described in [39] . For the present case we obtain 
the conversion formulae 


q( 1),finite 
^ ^Catani 
Q (2),finite 
^ ^Catani 


(-)(!),finite 
^(2),finite ^ 


(4.9) 


with All, in the case of a gg initial state, is given by 

All = -^7r^CA + i7r/3o- (4.10) 

In order to test the correctness of our results we have performed a number of checks, 
which we list in the following. 
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1. First of all, we computed explicitly all one- and two-loop diagrams relevant for 
gg —)• V 1 V 2 , including those diagrams in class B which are expected not to give 
any contribution due to Furry’s theorem, see Section 4. We have verified that, after 
reduction to master integrals, all diagrams in class B sum up to zero. 

2. We have verified explicitly that the coefficients Aj respect the expected symmetry 
relations derived in (2.27) and (2.28). 

3. We have verified explicitly that the IR poles of the two-loop amplitude have the 
structure predicted by Catani’s formula, see Section 4. This provides a strong check 
of the correctness of the result. 

4. We have performed a thorough comparison of our results with an independent cal¬ 

culation of the same process [58] . Specifically, we compared our results prior to UV 
renormalisation and IR subtraction. While the representation of the amplitudes in 
terms of spinor structures in [58] has a different form than our decomposition (3.12), 
we found that both are equivalent. For the full helicity amplitudes we have found 
perfect numerical agreement at one- and two-loop order. Moreover, expressing the 
form factors defined in [58] as linear combinations of our form factors we have 

verified that for each of them independently we have perfect numerical agreement at 
one- and two-loop order. 

5 Numerical C++ implementation and results 

For the numerical evaluation of the helicity amplitudes for gg —>■ V 1 V 2 —)• 4 leptons, we 
implemented our results for the form factors and at one- and two-loop order 

in a dedicated C++ code. The implementation is based on the solutions for the master inte¬ 
grals presented in [39], which were specifically constructed for fast and reliable numerical 
evaluations. We organised our form factor implementation in form of a library, which is 
supplemented by a simple command line interface. We provide the software package for 
public download on HepForge at http://vvamp.hepforge.org. 

For the numerical evaluation of the multiple polylogarithms encountered in the so¬ 
lutions for the master integrals, we employ their implementation [59] in the GiNaC [53] 
library. To identify and account for possible numerical instabilities of the form factors in 
collinear or other potentially problematic regions of phase space, the code compares nu¬ 
merical evaluations, which are obtained using different floating point data types, similar 
to the setup used in [39] . If the results obtained with different precision settings differ be¬ 
yond a user-defined tolerance, the code successively increases the precision until the target 
precision is met. 

For the rather central benchmark point of [40], the double precision mode of our code 
takes roughly 600ms on a single computer core and results in at least 11 significant digits 
for all of the Ej^^^ . In order to estimate the actual precision, the default behaviour of our 

code is to reevaluate the algebraic expressions in quad precision, which results in a total run¬ 
time of roughly 3s for this phase space point. The run-time can increase further for regions 
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Figure 2. Real parts of the two loop form factors for the process gg —>■ FlV 2 - The 

plots illustrate their dependence on the velocity, /Sa, and the cosine of the scattering angle, cos 03, 
of the vector boson Vi, where p\ = 2 p\ is chosen for the vector boson virtualities. 


close to the phase space boundaries, where the multiple polylogarithms take more time to 
evaluate and the precision control of our code may switch to higher precision computations 
in order to return reliable numbers. Depending on the precision setting and on the region of 
phase space, the evaluations of the multiple polylogarithms and of the algebraic coefficients 
may require comparable portions of the run-time. However, in double precision mode or 
close to the phase space boundaries, the run-time is dominated by multiple polylogarithm 
evaluations. The described version of our code implements a minimal set of 9 coefficients 
Aj and employs four evaluations of them with different kinematics in order to derive the 
remaining form factors using crossing relations. If required, it is straight-forward to further 
improve the evaluation speed, either by proper caching of multiple polylogarithms or, at 
the price of an increased code size, by an explicit implementation of all form factors, as we 
did for the process qq' —)• V 1 V 2 in [39]. 

In order to illustrate the form factors and the reliability of the code, we used the latter 
to plot the real part of the two-loop form factors for the case p\ = 2p| in Figures 2 and 3. In 
the plots, we vary the relativistic velocity (3^ and the cosine of the scattering angle cos ^3 of 
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Figure 3. Real parts of the two loop form factors for the process gg —>■ ViV2. The 

plots illustrate their dependence on the velocity, /Sa, and the cosine of the scattering angle, cos 03, 
of the vector boson Vi, where p\ = 2 p\ is chosen for the vector boson virtualities. 


the vector boson Vi, where /3s = K/{s+p‘^—p1) and cos 03 = (2t + s— p|)/k. Compared 
to the results for the form factors Ej in the process qq' —)• V 1 V 2 in [39], we observe strong 
enhancements for the forward, backward and production threshold regions for the form 
factors in the present case. However, for the physical helicity amplitudes (3.12) we wish 
to point out that an additional dampening (very) close to the aforementioned phase space 
boundaries should be taken into account due to the additional overall factors Cll and 
Clr (3.13). 


6 Conclusions 

In this paper we computed the two-loop massless QCD corrections to the helicity ampli¬ 
tudes for the production of pairs of off-shell electroweak gauge bosons, V 1 V 2 , in the gluon 
fusion channel. For the calculation we employed the solutions for the master integrals 
presented in [39]. Contracting the diboson amplitude with the leptonic decay currents we 
have constructed the helicity amplitudes for gg —)• V 1 V 2 —>• 4 leptons. We have compared 
our results to an independent calculation [58] and find perfect agreement. Our results for 
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these amplitudes provide the fundamental ingredient required to compute the NLO correc¬ 
tions to diboson production processes in gluon fusion. These corrections would contribute 
formally at N^LO to the processes pp ^ V 1 V 2 + X, but their inclusion may be important 
to match the expected experimental accuracy due to the large gluon luminosity at the 
LHC. In particular studying their impact is required to obtain a more reliable estimate of 
the theory uncertainty and to establish more precise constraints on the total Higgs decay 
width [20, 22, 23]. We provide both analytical results and a C++ code for the numerical 
evaluation of the amplitudes on HepForge at http://vvamp.hepforge.org. 
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A Form factor relations 


In this Appendix we present the explicit formulae needed in order to compute the 18 form 
factors defined for the amplitude (3.12), starting from the 20 form factors Aj defined 

in (2.13). For the Mllll amplitude we find 

cnLL _ + A2 + A3 Aiq 

-^1 , 2 2 ’ 
tu — S 

j^LL ^ ^ 14 (t - pI) - Ai2{u -pI)-S A4 ^ 

^ s{tu — pIp\) 2 s ’ 

pLL ^ ^ 14 (^ - Vl) - Ai3{u -pI)+A2+A3-S A5 ^ 

^ s{tu — p|p|) 2 s ’ 

t^LL _ — pI) — Ai2{t — P3) + A2+ A3 — sA^ Ai9 

^ s{tu — p|p|) 2 s ’ 


pLL ^ - pI) - Ai3{t -pI)-s A^ ^ 

^ s{tu — p\p\) 2 s ’ 


t^ll ^ ('» - A 3 ) Aio - Ai 4 


LL _ -P\){M - ^ 3 ) ^8 - Ai2 


= 


S{tu-plpi) 


t^ll ^ (•» -p\){M - A 3 ) Ag - Ai 3 
s situ-pjpl) + s 


E.LL {t - P3){^2 - A 3 ) All - Ai5 

-C/Q = -^-^- 


s{tu-plpl) 


(AT) 
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For the Mlrll amplitude we have instead 


pLR _ ^2+As ^16 

tu — P'1 S 




Et^ = - 


A 2 + A, 


Ai8 


= - 


s{tu-plpl) 

A 20 

’ 

{t - pI){A 2 + As) A 8 + A 12 


E^^ = - 




An 

Yi ’ 

A 2 + As 

s{tu-plpl) 


Aiq 

Yi ’ 


r^LR _ ~ ■P3)(^2 + ^ 3 ) _ ^10 + Ai4 

® s{tu — p‘^p‘l) S ’ 

{u - pI){A 2 + As) Ag + Ais 


s{tu-plpl) 


jpLR _ 

Es - - 


s{tu-plpl) 


LR _ {t — Ps){A2 + As) ^ 11 +^ 


TPEiJX _ 

Eg - 


15 


s{tu-pipi) 


(A.2) 


Let us consider the behaviour of the form factors Ej^^^ under the two permutations 7112 
and 7734 defined in (2.26). Using the crossing relations for the form factors A^^^ (2.27) and 
(2.28) one easily finds the corresponding ones for the form factors Ej^^^ To simplify our 
notation, we drop the superscript [A] in the following. Under permutation 7142 we obtain 

Ei^{s,u,pI,pI) = E^^{s,t,pl,pl), 

{u-pl)E^^{s,t,pl,pl) - {t-pl)E^^{s,t,pl,pl) 


E 2 {s,u,Ps,P4) = {s,t,Ps^P4:) + 

Ei^{s,u,pl,pl) = E^^{s,t,pl,pl) + 
Eq^{s,u,pI,pI) = -E^^{s,t,pl,pl) , 
and 

Ei^is,u,pl,pl) = E^^{s,t,pl,pl), 
Es^is,u,pl,pl) = E^^{s,t,pl,pl), 
E^^{s,u,pI,pI) = E^^{s,t,pl,pl). 


tu-plpl 


{t-pl)E^^{s,t,pl,pl) - {t-pl)E^^{s,t,pl,pl) 


tu-plpl 


E^^{s,u,pl,pi) = -EY{s,t,pl,pi) (A.3) 


'LL/ 


E^^{s,u,pI,pI) = E^^{s,t,pl,pl), 
E^^{s,u,pI,pI) = E^^{s,t,pl,pl), 


(A.4) 


Under permutation 7134 , instead, the form factors for both the LL and LR helicity config¬ 
urations transform in the same way 

u, pi pi) = pi, p1), 

Es''^^{s,U,pI,pI) = L^4^^is,t,pi,pi) , F;5l'^is,7i,pi,pi) = ^ 

776^1^^(5,11, pi,pi) = -E^^^^s,t,plpl), E^^^^{s,u,pIpI) = -S^^is,t,pi,pi). (A.5) 

Exploiting all of these crossing relations we find that only 9 out of the 18 form factors 
are effectively independent, while the other 9 can be obtained by the crossing rules above. 
The number of independent form factors Ej^^^ coincides with the number of independent 
form factors Aj found in Section 2. 
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